"""
一个类用来可视化OpenSeesPy模型
作者: 闫业祥
微信: yx592715024
目前博士就读 同济大学桥梁系 健康监测与振动控制教研室
"""

import numpy as np
import pyvista as pv
import openseespy.opensees as ops

# 如下为OpenSees中的单元类识别号，用来判断单元类型
ELE_TAG_Truss = (12, 13, 14, 15, 16, 17, 18, 169)  # 169 CatenaryCable
ELE_TAG_Link = (19, 20, 21, 22, 23, 24, 25, 26, 260, 27,  # zeroLength
                86,  # 86为 twoNodeLink
                84, 85, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104,
                105, 106, 107, 108, 109, 130, 131, 132, 147, 148, 149, 150, 151, 152, 153,
                158, 159, 160, 161, 165, 166)  # Bearing
ELE_TAG_Plane = (31, 32, 33, 40, 47, 50, 52, 53, 54, 55, 59, 61, 119, 120, 126, 134,
                 156, 157, 167, 168, 173, 174, 175, 203, 204, 207, 208, 209)
ELE_TAG_Joint = (71, 72, 81, 8181, 82, 83,)
ELE_TAG_Tetrahedron = (179,)
ELE_TAG_Brick = (36, 37, 38, 39, 41, 42, 43, 44, 45, 46, 48, 49, 51, 56, 57, 58,
                 121, 122, 127)


class model_base:
    """
    模型数据获取基类，不要实例化
    """

    def __init__(self):
        pass

    def _get_model_data(self):
        """
        用来获取有限元模型的数据
        :return: None
        """
        # Get all the node tags
        NodeTags = ops.getNodeTags()
        NodeTags.sort()
        num_node = len(NodeTags)
        # Get all the ele tags
        EleTags = ops.getEleTags()
        EleTags.sort()
        num_ele = len(EleTags)
        # Dict of node coordinates
        Node_Coords = np.zeros((num_node, 3))
        Node_Index = dict()  # key: nodeTag, value: index in Node_Coords
        model_dims = []

        for i, Tag in enumerate(NodeTags):
            coord = ops.nodeCoord(Tag)
            model_dim = len(coord)
            if model_dim == 1:
                coord.extend([0, 0])
            elif model_dim == 2:
                coord.extend([0])
            model_dims.append(model_dim)
            Node_Coords[i] = np.array(coord)
            Node_Index[Tag] = i
        points = Node_Coords

        # lines and faces
        # 单元的连接方式，节点数目n-节点号1-节点号2-...-节点号n
        truss_cells = []
        link_cells = []
        lines_cells = []
        all_lines_cells = []
        planes_cells = []
        tetrahedron_cells = []
        brick_cells = []
        ele_midpoints = []  # 单元中点坐标
        for i, ele in enumerate(EleTags):
            ele_nodes = ops.eleNodes(ele)
            # 以下根据单元节点数目来判断单元类型
            if len(ele_nodes) == 2:
                nodeI, nodeJ = ele_nodes
                idxI, idxJ = Node_Index[nodeI], Node_Index[nodeJ]
                all_lines_cells.append([2, idxI, idxJ])
                if ops.getEleClassTags(ele)[0] in ELE_TAG_Truss:
                    truss_cells.append([2, idxI, idxJ])
                elif ops.getEleClassTags(ele)[0] in ELE_TAG_Link:
                    link_cells.append([2, idxI, idxJ])
                else:
                    lines_cells.append([2, idxI, idxJ])
                ele_midpoints.append((Node_Coords[idxI] + Node_Coords[idxJ]) / 2)
            elif len(ele_nodes) == 3:
                nodeI, nodeJ, nodeK = ops.eleNodes(ele)
                idxI, idxJ, idxK = Node_Index[nodeI], Node_Index[nodeJ], Node_Index[nodeK]
                planes_cells.append([3, idxI, idxJ, idxK])
                ele_midpoints.append((Node_Coords[idxI] + Node_Coords[idxJ] + Node_Coords[idxK]) / 3)
            elif len(ele_nodes) == 4 or len(ele_nodes) == 9:
                if len(ele_nodes) == 4:
                    nodeI, nodeJ, nodeK, nodeL = ops.eleNodes(ele)
                else:
                    nodeI, nodeJ, nodeK, nodeL = ops.eleNodes(ele)[0:4]
                idxI, idxJ, idxK, idxL = Node_Index[nodeI], Node_Index[nodeJ], Node_Index[nodeK], Node_Index[nodeL]
                if ops.getEleClassTags(ele)[0] in ELE_TAG_Tetrahedron:  # 四面体
                    tetrahedron_cells.append([3, idxI, idxJ, idxK])
                    tetrahedron_cells.append([3, idxI, idxJ, idxL])
                    tetrahedron_cells.append([3, idxI, idxK, idxL])
                    tetrahedron_cells.append([3, idxJ, idxK, idxL])
                else:
                    planes_cells.append([4, idxI, idxJ, idxK, idxL])
                ele_midpoints.append(
                    (Node_Coords[idxI] + Node_Coords[idxJ] + Node_Coords[idxK] + Node_Coords[idxL]) / 4)
            elif len(ele_nodes) == 8 or len(ele_nodes) == 20:
                if len(ele_nodes) == 8:
                    node1, node2, node3, node4, node5, node6, node7, node8 = ops.eleNodes(ele)
                else:
                    node1, node2, node3, node4, node5, node6, node7, node8 = ops.eleNodes(ele)[0:8]
                tag_list = [Node_Index[node1], Node_Index[node2], Node_Index[node3], Node_Index[node4],
                            Node_Index[node5], Node_Index[node6], Node_Index[node7], Node_Index[node8]]
                temp_points = np.array([Node_Coords[i] for i in tag_list])
                idxxx = np.argsort(temp_points[:, 2])
                temp_points = temp_points[idxxx]
                tag_list = np.array(tag_list)[idxxx]
                temp_points = [tuple(i) for i in temp_points]
                tag_list = list(tag_list)
                tag_counter1 = counter_clockwise(temp_points[:4], tag_list[:4])  # 逆时针排序
                tag_counter2 = counter_clockwise(temp_points[4:], tag_list[4:])  # 逆时针排序
                tag_list = tag_counter1 + tag_counter2
                idx1, idx2, idx3, idx4, idx5, idx6, idx7, idx8 = tag_list
                brick_cells.append([4, idx1, idx4, idx3, idx2])
                brick_cells.append([4, idx5, idx6, idx7, idx8])
                brick_cells.append([4, idx1, idx5, idx8, idx4])
                brick_cells.append([4, idx2, idx3, idx7, idx6])
                brick_cells.append([4, idx1, idx2, idx6, idx5])
                brick_cells.append([4, idx3, idx4, idx8, idx7])
                ele_midpoints.append((Node_Coords[idx1] + Node_Coords[idx2] + Node_Coords[idx3] +
                                      Node_Coords[idx4] + Node_Coords[idx5] + Node_Coords[idx6] +
                                      Node_Coords[idx7] + Node_Coords[idx8]) / 8)
        ele_midpoints = np.array(ele_midpoints)

        # 自动根据模型节点坐标判断坐标轴边框的位置
        minNode = np.min(points, axis=0)
        maxNode = np.max(points, axis=0)
        space = (maxNode - minNode) / 15
        minNode = minNode - 2 * space
        maxNode = maxNode + 2 * space
        bounds = [minNode[0], maxNode[0], minNode[1], maxNode[1], minNode[2], maxNode[2]]
        max_bound = np.max(maxNode - minNode)

        # 有限元数据，包括点，线，面，实体的点连接
        self.model_data_finished = True
        self.points = points
        self.truss_cells = truss_cells
        self.link_cells = link_cells
        self.lines_cells = lines_cells
        self.all_lines_cells = all_lines_cells
        self.planes_cells = planes_cells
        self.tetrahedron_cells = tetrahedron_cells
        self.brick_cells = brick_cells
        self.ele_midpoints = ele_midpoints
        self.bounds = bounds
        self.max_bound = max_bound
        self.num_node = num_node
        self.num_ele = num_ele
        self.NodeTags = NodeTags
        self.EleTags = EleTags
        self.model_dims = model_dims

    def _eigen_analyze(self, modeTag, alpha):
        """
        本函数用来计算特征值，并保存相应数据
        :param modeTag: 标量，模态号
        :param alpha: 标量，缩放因子
        :return: None
        """
        mode_tag = modeTag
        ops.wipeAnalysis()
        if mode_tag == 1:
            eigenValues = ops.eigen(2)[:1]
        else:
            eigenValues = ops.eigen(mode_tag)
        omega = np.sqrt(eigenValues)
        f = omega / (2 * np.pi)

        if not self.model_data_finished:
            self._get_model_data()

        node_deform = {}
        eigen_vector = {}
        x_eigen = np.zeros(self.num_node)
        y_eigen = np.zeros(self.num_node)
        z_eigen = np.zeros(self.num_node)
        deform_points = np.zeros((self.num_node, 3))
        deform_anti_points = np.zeros((self.num_node, 3))

        if alpha is None:
            for i, Tag in enumerate(self.NodeTags):
                coord = ops.nodeCoord(Tag)
                eigen = ops.nodeEigenvector(Tag, mode_tag)
                if len(coord) == 1:
                    coord.extend([0, 0])
                    eigen.extend([0, 0])
                elif len(coord) == 2:
                    coord.extend([0])
                    eigen = eigen[:2]
                    eigen.extend([0])
                else:
                    eigen = eigen[:3]
                x_eigen[i], y_eigen[i], z_eigen[i] = eigen
            eigen_max = np.max(np.sqrt(x_eigen ** 2 + y_eigen ** 2 + z_eigen ** 2))
            alpha = self.max_bound / 10 / eigen_max

        for i, Tag in enumerate(self.NodeTags):
            coord = ops.nodeCoord(Tag)
            eigen = ops.nodeEigenvector(Tag, mode_tag)
            if len(coord) == 1:
                coord.extend([0, 0])
                eigen.extend([0, 0])
            elif len(coord) == 2:
                coord.extend([0])
                eigen = eigen[:2]
                eigen.extend([0])
            else:
                eigen = eigen[:3]
            x_eigen[i], y_eigen[i], z_eigen[i] = eigen
            node_deform[Tag] = np.array([eigen[ii] * alpha + coord[ii] for ii in range(3)])
            eigen_vector[Tag] = np.array(eigen)
            deform_points[i] = np.array([coord[ii] + eigen[ii] * alpha for ii in range(3)])
            deform_anti_points[i] = np.array([coord[ii] - eigen[ii] * alpha for ii in range(3)])

        self.mode_tag = mode_tag
        self.frequency = f
        self.eigen_points = deform_points
        self.eigen_anti_points = deform_anti_points
        self.xyz_eigen = np.sqrt(np.sum((self.eigen_points - self.points) ** 2, axis=1))


class OpenSeesVis(model_base):
    def __init__(self, point_size=6, line_width=6, colors_dict=None,
                 theme='paraview', notebook=False):
        """
        :param point_size: bool，是否显示节点编号
        :param line_width:  bool，是否显示单元编号
        :param colors_dict:  bool，是否显示坐标轴
        :param theme: str, 主题, ‘default', 'paraview', 'document', 'dark',
        :param notebook: bool, 是否嵌入在notebook中打开
        :return: None
        """
        super().__init__()
        self.point_size = point_size  # 节点尺寸
        self.line_width = line_width  # 线宽
        self.notebook = notebook
        # 设置主题
        pv.set_plot_theme(theme)

        self.model_data_finished = False  # 是否已提取模型数据
        self.mode_tag = None

        # 初始化颜色字典
        colors = dict(point='#840000', line='#0165fc', face='#06c2ac', solid='#f48924',
                      truss="#7552cc", link="#00c16e")
        if colors_dict is not None:
            colors.update(colors_dict)
        self.color_point = colors['point']
        self.color_line = colors['line']
        self.color_face = colors['face']
        self.color_solid = colors['solid']
        self.color_truss = colors['truss']
        self.color_link = colors['link']

    def model_vis(self, node_label=True, ele_label=True, outline=True, opacity=0.75):
        """
        本函数用来可视化模型，
        :param node_label: bool，是否显示节点编号
        :param ele_label:  bool，是否显示单元编号
        :param outline:  bool，是否显示坐标轴
        :param opacity:  float，面单元和体单元透明度
        :return: None
        """
        if not self.model_data_finished:
            self._get_model_data()

        plotter = pv.Plotter(notebook=self.notebook)
        # plotter.background_color = "white"

        point_plot = pv.PolyData(self.points)
        plotter.add_mesh(point_plot, color=self.color_point,
                         point_size=self.point_size, render_points_as_spheres=True)

        if len(self.truss_cells) > 0:
            truss_plot = pv.PolyData()
            truss_plot.points = self.points
            truss_plot.lines = np.array(self.truss_cells)
            plotter.add_mesh(truss_plot, color=self.color_truss,
                             render_lines_as_tubes=True, line_width=self.line_width)

        if len(self.link_cells) > 0:
            link_plot = pv.PolyData()
            link_plot.points = self.points
            link_plot.lines = np.array(self.link_cells)
            plotter.add_mesh(link_plot, color=self.color_link,
                             render_lines_as_tubes=False, line_width=self.line_width / 5)

        if len(self.lines_cells) > 0:
            line_plot = pv.PolyData()
            line_plot.points = self.points
            line_plot.lines = np.array(self.lines_cells)
            plotter.add_mesh(line_plot, color=self.color_line,
                             render_lines_as_tubes=True, line_width=self.line_width)

        if len(self.planes_cells) > 0:
            face_plot = pv.PolyData()
            face_plot.points = self.points
            face_plot.faces = np.hstack(self.planes_cells)
            plotter.add_mesh(face_plot, color=self.color_face, show_edges=True, opacity=opacity)

        if len(self.tetrahedron_cells) > 0:
            tet_plot = pv.PolyData()
            tet_plot.points = self.points
            tet_plot.faces = np.hstack(self.tetrahedron_cells)
            plotter.add_mesh(tet_plot, color=self.color_solid, show_edges=True, opacity=opacity)

        if len(self.brick_cells) > 0:
            bri_plot = pv.PolyData()
            bri_plot.points = self.points
            bri_plot.faces = np.hstack(self.brick_cells)
            plotter.add_mesh(bri_plot, color=self.color_solid, show_edges=True, opacity=opacity)

        plotter.add_text('OpenSees 3D View', position='upper_left',
                         font_size=15, color='black', font='courier', viewport=True)
        plotter.add_text('Num. of Node: {0} \n Num. of Ele:{1}'.format(self.num_node, self.num_ele),
                         position='upper_right', font_size=10, color='black', font='courier')
        if outline:
            plotter.show_bounds(grid=False, location='outer', bounds=self.bounds,
                                show_zaxis=True, color="black")
        if node_label:
            node_labels = [str(i) for i in self.NodeTags]
            plotter.add_point_labels(self.points, node_labels, text_color='white',
                                     font_size=12, bold=False, always_visible=True)
        if ele_label:
            ele_labels = [str(i) for i in self.EleTags]
            plotter.add_point_labels(self.ele_midpoints, ele_labels, text_color='#ff796c',
                                     font_size=12, bold=False, always_visible=True)
        plotter.add_axes(color='black')
        plotter.view_isometric()
        if np.max(self.model_dims) <= 2:
            plotter.view_xy(negative=False)
        plotter.show()

    def eigen_vis(self, modeTag=1, alpha=None, outline=False,
                  show_origin=False, colormap='coolwarm', opacity=1):
        """
        :param modeTag: int, 模态号
        :param alpha: float, 缩放因子，默认值根据最大变形取模型边界的1/5
        :param outline: bool, 是否显示坐标轴
        :param show_origin: bool, 是否显示原始形状
        :param colormap: 云图类型, 常用有 jet, rainbow, hot, afmhot, copper,
                            winter, cool, coolwarm, gist_earth, bone, binary, gray
        :param opacity: float, 面和体单元的透明度
        """

        self._eigen_analyze(modeTag, alpha)

        plotter = pv.Plotter(notebook=self.notebook)

        point_plot = pv.PolyData(self.eigen_points)
        plotter.add_mesh(point_plot, colormap=colormap, scalars=self.xyz_eigen,
                         interpolate_before_map=True, point_size=self.point_size,
                         render_points_as_spheres=True, smooth_shading=True)

        if len(self.all_lines_cells) > 0:
            if show_origin:
                line_plot_origin = pv.PolyData()
                line_plot_origin.points = self.points
                line_plot_origin.lines = np.array(self.all_lines_cells)
                plotter.add_mesh(line_plot_origin, color='gray',
                                 render_lines_as_tubes=False, line_width=self.line_width / 3)
            line_plot = pv.PolyData()
            line_plot.points = self.eigen_points
            line_plot.lines = np.array(self.all_lines_cells)
            plotter.add_mesh(line_plot, colormap=colormap, scalars=self.xyz_eigen,
                             interpolate_before_map=True, smooth_shading=True,
                             render_lines_as_tubes=True, line_width=self.line_width)
        if len(self.planes_cells) > 0:
            if show_origin:
                face_plot_origin = pv.PolyData()
                face_plot_origin.points = self.points
                face_plot_origin.faces = np.hstack(self.planes_cells)
                plotter.add_mesh(face_plot_origin, color='gray', style='wireframe',
                                 show_edges=True, line_width=self.line_width / 3)
            face_plot = pv.PolyData()
            face_plot.points = self.eigen_points
            face_plot.faces = np.hstack(self.planes_cells)
            plotter.add_mesh(face_plot, colormap=colormap, scalars=self.xyz_eigen,
                             show_edges=True, opacity=opacity, smooth_shading=True,
                             show_scalar_bar=True, interpolate_before_map=True)
        if len(self.tetrahedron_cells) > 0:
            if show_origin:
                tet_plot_origin = pv.PolyData()
                tet_plot_origin.points = self.points
                tet_plot_origin.faces = np.hstack(self.tetrahedron_cells)
                plotter.add_mesh(tet_plot_origin, color='gray', style='wireframe',
                                 show_edges=True, line_width=self.line_width / 3)
            tet_plot = pv.PolyData()
            tet_plot.points = self.eigen_points
            tet_plot.faces = np.hstack(self.tetrahedron_cells)
            plotter.add_mesh(tet_plot, colormap=colormap, scalars=self.xyz_eigen,
                             show_edges=True, opacity=opacity, smooth_shading=True,
                             show_scalar_bar=True, interpolate_before_map=True)

        if len(self.brick_cells) > 0:
            if show_origin:
                bri_plot_origin = pv.PolyData()
                bri_plot_origin.points = self.points
                bri_plot_origin.faces = np.hstack(self.brick_cells)
                plotter.add_mesh(bri_plot_origin, color='gray', style='wireframe',
                                 show_edges=True, line_width=self.line_width / 3)
            bri_plot = pv.PolyData()
            bri_plot.points = self.eigen_points
            bri_plot.faces = np.hstack(self.brick_cells)
            plotter.add_mesh(bri_plot, colormap=colormap, scalars=self.xyz_eigen,
                             show_edges=True, opacity=opacity, smooth_shading=True,
                             show_scalar_bar=True, interpolate_before_map=True)

        plotter.add_scalar_bar(color='#000000', n_labels=10, label_font_size=12)

        plotter.add_text('OpenSees 3D View', position='upper_left',
                         font_size=16, color='black', font='courier')
        plotter.add_text('Mode {}\nf = {:.3f} Hz\nT = {:.3f} s'.format(self.mode_tag,
                                                                       self.frequency[-1], 1 / self.frequency[-1]),
                         position='upper_right', font_size=12, color='black', font='courier')
        if outline:
            plotter.show_bounds(grid=False, location='outer', bounds=self.bounds,
                                show_zaxis=True, color="black")
        plotter.add_axes(color='black')
        plotter.view_isometric()
        if np.max(self.model_dims) <= 2:
            plotter.view_xy(negative=False)
        plotter.show()

    def eigen_animation(self, output_file, modeTag=1, alpha=None, outline=False,
                        colormap='coolwarm', opacity=1):
        """
        :param output_file: str, 输出文件字符串
        :param modeTag: int, 模态号
        :param alpha: float, 缩放因子，默认值根据最大变形取模型边界的1/5
        :param outline: bool, 是否显示坐标轴
        :param colormap: 云图类型, 常用有 jet, rainbow, hot, afmhot, copper,
                            winter, cool, coolwarm, gist_earth, bone, binary, gray
        :param opacity: float, 面和体单元的透明度
        """

        self._eigen_analyze(modeTag, alpha)

        plotter = pv.Plotter(notebook=self.notebook)

        point_plot = pv.PolyData(self.points)
        plotter.add_mesh(point_plot, colormap=colormap, scalars=self.xyz_eigen,
                         point_size=self.point_size, render_points_as_spheres=True,
                         interpolate_before_map=True)

        if len(self.all_lines_cells) > 0:
            line_plot = pv.PolyData()
            line_plot.points = self.points
            line_plot.lines = np.array(self.all_lines_cells)
            plotter.add_mesh(line_plot, colormap=colormap, scalars=self.xyz_eigen,
                             render_lines_as_tubes=True, line_width=self.line_width,
                             interpolate_before_map=True)
        if len(self.planes_cells) > 0:
            face_plot = pv.PolyData()
            face_plot.points = self.points
            face_plot.faces = np.hstack(self.planes_cells)
            plotter.add_mesh(face_plot, colormap=colormap, scalars=self.xyz_eigen,
                             show_edges=True, opacity=opacity,
                             show_scalar_bar=True, interpolate_before_map=True)
        if len(self.tetrahedron_cells) > 0:
            tet_plot = pv.PolyData()
            tet_plot.points = self.points
            tet_plot.faces = np.hstack(self.tetrahedron_cells)
            plotter.add_mesh(tet_plot, colormap=colormap, scalars=self.xyz_eigen,
                             show_edges=True, opacity=opacity,
                             show_scalar_bar=True, interpolate_before_map=True)

        if len(self.brick_cells) > 0:
            bri_plot = pv.PolyData()
            bri_plot.points = self.points
            bri_plot.faces = np.hstack(self.brick_cells)
            plotter.add_mesh(bri_plot, colormap=colormap, scalars=self.xyz_eigen,
                             show_edges=True, opacity=opacity,
                             show_scalar_bar=True, interpolate_before_map=True)

        plotter.add_text('OpenSees 3D View', position='upper_left',
                         font_size=16, color='black', font='courier')
        plotter.add_text('Mode {}\nf = {:.3f} Hz\nT = {:.3f} s'.format(self.mode_tag,
                                                                       self.frequency[-1], 1 / self.frequency[-1]),
                         position='upper_right', font_size=12, color='black', font='courier')
        plotter.add_scalar_bar(color='#000000', n_labels=10, label_font_size=12)
        if outline:
            plotter.show_bounds(grid=False, location='outer', bounds=self.bounds,
                                show_zaxis=True, color="black")
        plotter.add_axes(color='black')
        plotter.view_isometric()
        if np.max(self.model_dims) <= 2:
            plotter.view_xy(negative=False)

        # 动画
        plotter.open_gif(output_file)
        plt_points = [self.eigen_anti_points, self.points, self.eigen_points]
        index = [1, 2, 1, 0] * 10
        for idx in index:
            points = plt_points[idx]
            xyz = self.points - points
            xyz_eigen = np.sqrt(np.sum(xyz ** 2, axis=1))
            plotter.update_coordinates(points, mesh=point_plot, render=False)
            if len(self.all_lines_cells) > 0:
                plotter.update_coordinates(points, mesh=line_plot, render=False)
                plotter.update_scalars(xyz_eigen, mesh=line_plot, render=False)
            if len(self.planes_cells) > 0:
                plotter.update_coordinates(points, mesh=face_plot, render=False)
                plotter.update_scalars(xyz_eigen, mesh=face_plot, render=False)
            if len(self.tetrahedron_cells) > 0:
                plotter.update_coordinates(points, mesh=tet_plot, render=False)
                plotter.update_scalars(xyz_eigen, mesh=tet_plot, render=False)
            if len(self.brick_cells) > 0:
                plotter.update_coordinates(points, mesh=bri_plot, render=False)
                plotter.update_scalars(xyz_eigen, mesh=bri_plot, render=False)
            plotter.update_scalar_bar_range(clim=[np.min(xyz_eigen), np.max(xyz_eigen)], name=None)
            plotter.write_frame()
        plotter.show()
        plotter.close()


def sort_xy(x, y):
    # 逆时针排序点
    x0 = np.mean(x)
    y0 = np.mean(y)
    r = np.sqrt((x - x0) ** 2 + (y - y0) ** 2)
    angles = np.where((y - y0) >= 0, np.arccos((x - x0) / r),
                      4 * np.pi - np.arccos((x - x0) / r))
    mask = np.argsort(angles)
    x_max = np.max(x)
    if x[mask][0] != x_max:
        mask = np.roll(mask, 1)
    return mask


def counter_clockwise(points, tag):
    """
    用来将一个面上的点逆时针排序
    """
    x = np.array([point[0] for point in points])
    y = np.array([point[1] for point in points])
    z = np.array([point[2] for point in points])

    if all(np.abs(x - x[0]) < 1e-6):  # yz平面
        # def algo(point):
        #    return (math.atan2(point[2] - z_center, point[1] - y_center) + 2 * math.pi) % (2*math.pi)
        # sorted_points = sorted(points,key = algo )
        index = sort_xy(y, z)
        sorted_points = list(zip(x[index], y[index], z[index]))
        sorted_id = [points.index(i) for i in sorted_points]
        sorted_tag = [tag[id] for id in sorted_id]
    elif all(np.abs(y - y[0]) < 1e-6):  # xz平面
        # def algo(point):
        #    return (math.atan2(point[2] - z_center, point[0] - x_center) + 2 * math.pi) % (2*math.pi)
        # sorted_points = sorted(points,key = algo )
        index = sort_xy(x, z)
        sorted_points = list(zip(x[index], y[index], z[index]))
        sorted_id = [points.index(i) for i in sorted_points]
        sorted_tag = [tag[id] for id in sorted_id]
    else:
        # def algo(point):
        #    return (math.atan2(point[1] - y_center, point[0] - x_center) + 2 * math.pi) % (2*math.pi)
        # sorted_points = sorted(points,key = algo )
        index = sort_xy(x, y)
        sorted_points = list(zip(x[index], y[index], z[index]))
        sorted_id = [points.index(i) for i in sorted_points]
        sorted_tag = [tag[i] for i in sorted_id]
    return sorted_tag


if __name__ == "__main__":

    def model():
        pi = 3.1415

        AnalysisType = "Pushover"  # Cyclic   Pushover   Gravity

        ops.wipe()

        ops.model('basic', '-ndm', 3, '-ndf', 6)

        ###################################
        # Define Material
        ###################################

        # Define PSUMAT and convert it to plane stress material
        ops.nDMaterial('PlaneStressUserMaterial', 1, 40, 7, 20.7e6, 2.07e6, -4.14e6, -0.002, -0.01, 0.001, 0.3)
        ops.nDMaterial('PlateFromPlaneStress', 4, 1, 1.25e10)

        # Define material for rebar
        ops.uniaxialMaterial('Steel02', 7, 379e6, 202.7e9, 0.01, 18.5, 0.925, 0.15)
        ops.uniaxialMaterial('Steel02', 8, 392e6, 200.6e9, 0.01, 18.5, 0.925, 0.15)

        # Convert rebar material to plane stress/plate rebar
        # Angle 0 is for vertical rebar and 90 is for horizontal rebar
        ops.nDMaterial('PlateRebar', 9, 7, 90.0)
        ops.nDMaterial('PlateRebar', 10, 8, 90.0)
        ops.nDMaterial('PlateRebar', 11, 8, 0.0)

        # Define LayeredShell sections. Section 1 is used for the special boundary elements and section 2 is used for
        # the unconfined interior wall portion
        ops.section('LayeredShell', 1, 10, 4, 0.0125, 11, 0.0002403, 11, 0.0003676, 4, 0.024696, 4, 0.024696, 4,
                    0.024696, 4,
                    0.024696, 11, 0.0003676, 11, 0.0002403, 4, 0.0125)
        ops.section('LayeredShell', 2, 8, 4, 0.0125, 11, 0.0002403, 10, 0.0002356, 4, 0.0495241, 4, 0.0495241, 10,
                    0.0002356,
                    11, 0.0002403, 4, 0.0125)

        # ##################
        # NODES
        # ##################
        # define nodes
        for i in range(11):
            ops.node(16 * i + 1, 0, 0, 0.2 * i)
            ops.node(16 * i + 2, 0.2, 0, 0.2 * i)
            ops.node(16 * i + 3, 0.5, 0, 0.2 * i)
            ops.node(16 * i + 4, 0.8, 0, 0.2 * i)
            ops.node(16 * i + 5, 1.0, 0, 0.2 * i)

            ops.node(16 * i + 6, 0, 1.0, 0.2 * i)
            ops.node(16 * i + 7, 0.2, 1, 0.2 * i)
            ops.node(16 * i + 8, 0.5, 1, 0.2 * i)
            ops.node(16 * i + 9, 0.8, 1, 0.2 * i)
            ops.node(16 * i + 10, 1.0, 1, 0.2 * i)

            ops.node(16 * i + 11, 0, 0.2, 0.2 * i)
            ops.node(16 * i + 12, 0, 0.5, 0.2 * i)
            ops.node(16 * i + 13, 0, 0.8, 0.2 * i)
            ops.node(16 * i + 14, 1, 0.2, 0.2 * i)
            ops.node(16 * i + 15, 1, 0.5, 0.2 * i)
            ops.node(16 * i + 16, 1, 0.8, 0.2 * i)

        ##########################
        # ELEMENTS
        ##########################

        ShellType = "ShellNLDKGQ"
        # ShellType = "ShellMITC4"
        for i in range(10):
            ops.element(ShellType, i + 1, 16 * i + 1, 16 * i + 2, 16 * (i + 1) + 2, 16 * (i + 1) + 1, 1)
            ops.element(ShellType, i + 11, 16 * i + 2, 16 * i + 3, 16 * (i + 1) + 3, 16 * (i + 1) + 2, 2)
            ops.element(ShellType, i + 21, 16 * i + 3, 16 * i + 4, 16 * (i + 1) + 4, 16 * (i + 1) + 3, 2)
            ops.element(ShellType, i + 31, 16 * i + 4, 16 * i + 5, 16 * (i + 1) + 5, 16 * (i + 1) + 4, 1)

            ops.element(ShellType, i + 41, 16 * i + 6, 16 * i + 7, 16 * (i + 1) + 7, 16 * (i + 1) + 6, 1)
            ops.element(ShellType, i + 51, 16 * i + 7, 16 * i + 8, 16 * (i + 1) + 8, 16 * (i + 1) + 7, 2)
            ops.element(ShellType, i + 61, 16 * i + 8, 16 * i + 9, 16 * (i + 1) + 9, 16 * (i + 1) + 8, 2)
            ops.element(ShellType, i + 71, 16 * i + 9, 16 * i + 10, 16 * (i + 1) + 10, 16 * (i + 1) + 9, 1)

            ops.element(ShellType, i + 81, 16 * i + 1, 16 * i + 17, 16 * i + 27, 16 * i + 11, 1)
            ops.element(ShellType, i + 91, 16 * i + 11, 16 * i + 27, 16 * i + 28, 16 * i + 12, 2)
            ops.element(ShellType, i + 101, 16 * i + 12, 16 * i + 28, 16 * i + 29, 16 * i + 13, 2)
            ops.element(ShellType, i + 111, 16 * i + 13, 16 * i + 29, 16 * i + 22, 16 * i + 6, 1)

            ops.element(ShellType, i + 121, 16 * i + 5, 16 * i + 14, 16 * i + 30, 16 * i + 21, 1)
            ops.element(ShellType, i + 131, 16 * i + 14, 16 * i + 15, 16 * i + 31, 16 * i + 30, 2)
            ops.element(ShellType, i + 141, 16 * i + 15, 16 * i + 16, 16 * i + 32, 16 * i + 31, 2)
            ops.element(ShellType, i + 151, 16 * i + 16, 16 * i + 10, 16 * i + 26, 16 * i + 32, 1)

        # P-delta columns

        for i in range(10):
            ops.element('truss', 500 + i + 1, 16 * i + 1, 16 * (i + 1) + 1, 223.53e-6, 7)
            ops.element('truss', 500 + i + 11, 16 * i + 2, 16 * (i + 1) + 2, 223.53e-6, 7)
            ops.element('truss', 500 + i + 21, 16 * i + 4, 16 * (i + 1) + 4, 223.53e-6, 7)
            ops.element('truss', 500 + i + 31, 16 * i + 5, 16 * (i + 1) + 5, 223.53e-6, 7)

            ops.element('truss', 500 + i + 41, 16 * i + 6, 16 * (i + 1) + 6, 223.53e-6, 7)
            ops.element('truss', 500 + i + 51, 16 * i + 7, 16 * (i + 1) + 7, 223.53e-6, 7)
            ops.element('truss', 500 + i + 61, 16 * i + 9, 16 * (i + 1) + 9, 223.53e-6, 7)
            ops.element('truss', 500 + i + 71, 16 * i + 10, 16 * (i + 1) + 10, 223.53e-6, 7)

            ops.element('truss', 500 + i + 81, 16 * i + 11, 16 * (i + 1) + 11, 223.53e-6, 7)
            ops.element('truss', 500 + i + 91, 16 * i + 13, 16 * (i + 1) + 13, 223.53e-6, 7)

            ops.element('truss', 500 + i + 101, 16 * i + 14, 16 * (i + 1) + 14, 223.53e-6, 7)
            ops.element('truss', 500 + i + 111, 16 * i + 16, 16 * (i + 1) + 16, 223.53e-6, 7)

        for i in range(176):
            m = 246000 * 4 / 9.8 / 176
            ops.mass(i + 1, m, m, m, 0, 0, 0)
        # Fix all bottom nodes
        ops.fixZ(0, 1, 1, 1, 1, 1, 1)


    model()

    A = OpenSeesVis(theme='paraview')
    A.model_vis()
    A.eigen_vis(modeTag=1)
    A.eigen_animation(output_file=r'C:\Users\yan\Desktop\yan.gif', modeTag=1, alpha=None, outline=False,
                      colormap='coolwarm', opacity=1)
    # A.static_vis(num_steps=10, alpha=500, node_label=False, ele_label=False, show_origin=False)
    # A.history_gif(r'C:\Users\yan\Desktop\palte.mp4', num_steps=400, dt=0.1, alpha=100)
